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ABSTRACT 

MicroRNAs (miRNAs) are ubiquitously expressed 
small non-coding RNAs that, in most cases, 
negatively regulate gene expression at the post- 
transcriptional level. miRNAs are involved in 
fine-tuning fundamental cellular processes such as 
proliferation, cell death and cell cycle control and 
are believed to confer robustness to biological 
responses. Here, we investigated simultaneously 
the transcriptional changes of miRNA and mRNA 
expression levels over time after activation of the 
Janus kinase/Signal transducer and activator of 
transcription (Jak/STAT) pathway by interferon-y 
stimulation of melanoma cells. To examine global 
miRNA and mRNA expression patterns, time-series 
microarray data were analysed. We observed 
delayed responses of miRNAs (after 24-48 h) with 
respect to mRNAs (12-24 h) and identified biological 
functions involved at each step of the cellular 
response. Inference of the upstream regulators 
allowed for identification of transcriptional regula- 
tors involved in cellular reactions to interferon-y 
stimulation. Linking expression profiles of transcrip- 
tional regulators and miRNAs with their annotated 
functions, we demonstrate the dynamic interplay of 
miRNAs and upstream regulators with biological 
functions. Finally, our data revealed network 
motifs in the form of feed-forward loops involving 
transcriptional regulators, mRNAs and miRNAs. 
Additional information obtained from integrating 
time-series mRNA and miRNA data may represent 
an important step towards understanding the 
regulatory principles of gene expression. 



INTRODUCTION 

MicroRNAs (miRNAs) have been discovered in 1993, and 
initially, these small non-coding RNAs have not attracted 
much interest from the scientific community (1). However, 
in recent years, it has emerged that the highly conserved 
and ubiquitously expressed miRNAs are of paramount 
importance for the regulation of gene expression in 
humans, animals and plants (2). Thus far, >1600 mature 
miRNAs have been identified in humans (mirBase version 
19), and each miRNA is predicted to regulate several 
hundreds of target genes, leading to the conservative 
estimate of >60% of human protein-coding genes being 
regulated by miRNAs (3,4). The binding of miRNAs to 
their target mRNAs generally results in mRNA down- 
regulation or degradation, with subsequent repression of 
protein synthesis (2,5). A common and estabhshed feature 
is that miRNAs do not need an entirely complementary 
region in the 3^ UTR of the target gene mRNA to bind to 
but can do with varying numbers of mismatching nucleo- 
tides. This renders in silico predictions of miRNA target 
genes very difficult, and thus far no efficient algorithm 
exists, which is able to reliably predict all, but no false- 
positive, target genes (6). 

Given the large number of protein-encoding genes 
that miRNAs can regulate post-transcriptionally, it is 
evident that they modulate and fine-tune almost all biolo- 
gical processes (7). Consequently, miRNAs have been 
impHcated in the regulation of processes that promote 
cancer growth, or conversely, in processes that might 
prevent cancers and other diseases from developing 
(8-11). Considering their tremendous regulatory potential 
and their often tissue- and disease-specific expression 
patterns (12), de-regulated individual miRNAs or altered 
global miRNA expression profiles could be indicative of 
disease risks and burdens; therefore, miRNAs are currently 
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being assessed as possible biomarkers to aid diagnosis 
and prediction of different types and stages of cancers, 
including melanoma (13,14). In addition, miRNAs are dis- 
cussed as targets for cancer therapeutics and as possible 
biomarkers (15). Despite recent progress in understanding 
miRNA effects on cell behaviour, the precise mechanisms 
and impHcations of miRNA actions are currently debated. 
To answer these questions, the dynamic regulation of 
miRNA expression changes will have to be considered, 
which thus far has been largely neglected (16). 

The initial point of regulation of miRNA biogenesis, the 
transcription of miRNA genes, is a tightly controlled 
multi-step process, which often involves auto-regulatory 
feedback loops and feed-forward loops (FFLs) in which 
miRNAs participate together with transcription factors 
(TFs) (7,17-19). Gene expression, in general, results in 
variable levels of gene transcripts and proteins. Together 
with expression noise, the magnitude of which is influenced 
by intrinsic and extrinsic factors (20), gene transcription 
and its inferred regulatory networks can be considered as 
'dynamic information processing systems' (21-23). 
However, the fluctuations in transcript levels (expression 
noise) have to be counter-balanced by a certain level of 
robustness in the biological responses, and this sturdiness 
is thought to be maintained by miRNAs (24). In this 
context, the integration of matching mRNA and miRNA 
data sets will become increasingly important. Recently, 
Muniategui et al. (25) have reviewed and grouped math- 
ematical and computational approaches for analysing the 
interplay between miRNAs and mRNA into three main 
categories: dependency analysis, linear regression and 
Bayesian methods. It was further emphasized that models 
combining heterogeneous experimental data, such as 
time-series data, would be more reliable to predict 
miRNA-mRNA interactions. Dynamic data of a given 
biological system can add valuable information to a 
better understanding of the underlying cellular processes 
that might be missed using cross-sectional data that only 
focus on single time points (26). Recently, Kim et al. (27) 
have analysed complex network dynamics by using 
time-series-derived expression data; with principal 
network analysis, they were able to capture major activa- 
tion patterns from two data sets and to generate the 
associated sub-networks and their interactions. Jayaswal 
et al. (28), in contrast, used odds ratio statistics to 
perform an integrative analysis on matched miRNA and 
mRNA time-course microarray data, which identified 
miRNAs with regulatory potential and their targeted 
mRNAs. Associations between TFs and miRNAs in 
monocytic differentiation were also determined in a time- 
lagged expression correlation analysis, which identified 
12 TFs regulating the expression of several key miRNAs 
(29). The importance of time-series gene expression data 
was also underscored by a recent review in which 
Bar- Joseph et al. (30) expertly summarized current know- 
ledge on this topic: different biological scenarios lead 
to different response patterns or programs, resulting 
in cyclic, sustained or most commonly peaked impulse 
responses after a stimulus and/or environmental 
perturbations. 



To investigate whether integrative time- series-derived 
data would provide a means to better explain and 
identify complex regulatory interactions, we generated 
data sets representing a melanoma cell-derived 
miRNome and transcriptome (mRNAs) analysed at dif- 
ferent time points after a transcriptional activation 
stimulus. We developed an analysis pipeline and 
combined known methods to extract information from 
these dynamic data sets, aiming at the visualization of 
functional variations that are connected to expression 
changes. We stimulated melanoma cells with interferon-y 
(IFN-y), a type II cytokine, which is known to induce 
STAT 1 -mediated growth inhibition and anti-proliferative 
effects in these cells (31,32). We set out to find potential 
explanations for these biological effects by integration of 
dynamic miRNA and mRNA data sets. Time-series dif- 
ferential expression analyses were performed, mainly in 
the form of contrasts between experimental conditions 
(IFN-y-treated samples versus untreated controls) using 
the R/Bioconductor package 'limma' (33), in combination 
with profile correlation analysis. Ingenuity® Pathway 
Analysis (IPA) and data visualization with Circos (34). 
We and others have shown before that activation of 
STAT TFs specifically up-regulates the expression of 
several miRNAs in a time- and dose-dependent manner 
(35-40). In the current study, we have observed a delayed 
response of the miRNome with respect to the transcrip- 
tome, and have dynamically connected biological func- 
tions involved in cellular responses to IFN-y. The 
up-regulation of several TFs downstream of STATl, 
which might be required for full transcriptional activation 
of most miRNAs as well as different processing times of 
primary miRNA transcripts (41), might account for this 
delayed miRNome response. In addition, we have 
identified several FFLs including TFs, miRNAs and 
mRNAs. We argue that the integration of time-series- 
derived miRNA and mRNA expression data provides 
valuable information for generation of biological 
networks and is highly relevant for fully understanding 
the regulation of biological responses and processes. 

MATERIALS AND METHODS 

Time-series experiments 

To investigate the dynamic regulation and potential 
co-regulation of miRNAs and mRNAs, time-course 
experiments were performed as described before (40). 
An overview of the experimental setup is shown in 
Supplementary Figure SI. Briefly, the time-course experi- 
ments were carried out using the human A375 melanoma 
cell line (ATCC, CRL-1619™). Cells were seeded into six- 
well plates at a density of 50 000 cells/well and were either 
left untreated (ctrl) or stimulated with human IFN-y 
(PeproTech Inc., final concentration of 50ng/mL) for 
the indicated time points. In parallel, a second control 
time-series experiment was performed including a 
pre-treatment step with 5 |iM Janus kinase (Jak) inhibitor 
I (JII; Calbiochem), added 1 h before commencing IFN-y 
stimulation. For microarray applications, RNA was 
extracted using the miRNeasy Mini kit (Qiagen); for 



Nucleic Acids Research, 2013, Vol. 41, No. 5 2819 



quantitative real-time PCR (qPCR) validations, TRIsure 
(Bioline) was used. The description of the qPCR method 
and primer sequences is given in Supplementary Data 
(Supplementary Table SI). We extracted and analysed 
RNAs from three biological repHcates, each consisting 
of three technical replicates. RNA quality was assessed 
using Nanodrop 2000 Spectrophotometer (Thermo 
Scientific) and Agilent 2100 Bioanalyzer (Agilent 
Technologies). 

miRNA and mRNA microarray expression profiling 

miRNA microarray s were performed using the Affymetrix 
GeneChip® miRNA 2.0 Array technology as described 
by Reinsbach et al. (40). In addition, matching mRNA 
microarrays for selected time points (ctrl, 3, 12, 24, 48 
and 72 h) were performed using Gene Chip® Human 
Gene 1.0 ST Array (Affymetrix, Santa Clara, CA, USA) 
on 17 samples. miRNA and mRNA microarray expres- 
sion data are available at ArrayExpress (www.ebi.ac.uk/ 
array express) under accession numbers E-MEXP-3544 
and E-MEXP-3720, respectively. 

Microarray data analysis 

The workflow of data processing and analysis is outlined 
in Figure 1. Pre-processing of microarray data was per- 
formed with Partek® Genomics Suite version 6.5 (Partek® 



GS) using the robust multi-chip analysis algorithm, which 
performs background adjustment, quantile normalization 
and probe summarization as described before (42). 
For mRNA data, GC content correction was used, as sug- 
gested by the default pipeHne of Partek® GS. To decrease 
the number of uninformative features, we used a filtering 
step: features, for which the maximum expression over all 
arrays did not reach a signal intensity of 7 in a log2 scale, 
were removed. Different methods ['betr' (43), 'timecourse' 
(44) and 'limma' (33)] to identify significantly differentially 
expressed (SDE) miRNAs and mRNAs were tested 
(Supplementary Data and Supplementary Table S2). 
In this study, linear models and the empirical Bayes 
method from the 'limma' package of R/Bioconductor 
were used for differential analysis of the miRNA and 
mRNA time-series data sets as described before (45). In 
addition, to obtain a Hst of features significantly regulated 
at each time point, we used the same 'limma' model and a 
set of contrasts, comparing expression in IFN-y-treated 
samples versus untreated controls. To control for false 
discovery rates (FDR), we used the Benjamini-Hochberg 
method to adjust P-values. 

Microarray data were visualized as heat maps using 
a standard eponymous function of R. A hierarchical 
clustering algorithm was applied to determine similar 
patterns in miRNA and mRNA expression profiles. 
Clustering was supported by bootstrapping using the 
'clusterCons' package of R (46). 
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Figure 1. Graphical representation of the computational workflow. In 
the first step, miRNA and mRNA array data were pre-processed 
and the quahty was assessed using Partek® GS. A filtering step was 
included to select only expressed and detectable features. Co-expression 
analysis and identification of potential targets were performed directly 
on the paired miRNA-mRNA data. Differential expression analysis 
was performed using the 'limma' package of R/Bioconductor to 
identify significant differentially expressed mRNAs and miRNAs over 
time. Functional and pathway analysis together with identification of 
upstream regulators was carried out in IPA based on SDE molecules. 
The dynamic behaviour of gene regulation circuitry in response to 
IFN-y stimulus was visualized as Circos plots. 



Correlation analyses using CoExpress 

Negatively correlated miRNA-mRNA pairs were deter- 
mined by the in-house developed stand-alone software 
tool CoExpress (freely available at www.bioinformatics. 
lu/CoExpress), which performs interactive detection 
of correlated profiles in large expression data sets. 
The software allows single-data type analysis for identifi- 
cation of miRNA-mRNA co-expression events as well as 
two-data type analysis for detection of correlated miRNA- 
mRNA expressions. Details on CoExpress are given in the 
Supplementary Data section. 

Data mining and visualization 

Canonical pathway analysis. Lists of mRNAs and 
miRNAs, differentially expressed between each condition 
(time points) and untreated controls (with FDR < 0.001), 
were uploaded in the IPA tool (Ingenuity® Systems, www. 
ingenuity.com) and analysed based on the IPA library of 
canonical pathways (content date 2012-05-08). The signifi- 
cance of the association between each list and a canonical 
pathway was measured by Fisher's exact test. As a result, 
a P-value was obtained, determining the probability that 
the association between the genes in our data set and a 
canonical pathway can be explained by chance alone. 
Functional analysis. To identify biological functions and/ 
or diseases that were most significant to our data sets 
functional analysis was done. As in canonical pathway 
analysis, features that met the FDR cut-off of 0.001 
when comparing each condition with the untreated 
control were analysed by IPA. Right-tailed Fisher's 
exact test was used to calculate a significant P-value for 
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each functional category as referenced in Ingenuity® 
Knowledge Base. The obtained P- value was further 
adjusted using the Benjamini-Hochberg correction. We 
focused our analysis on IPA categories with an adjusted 
P<0.01 (presented in at least one time point) directly 
related to cellular functions and diseases, providing dy- 
namical information about enriched biological functions. 
Upstream regulator analysis. Based on previous knowledge 
of expected effects between transcriptional regulators 
(TRs) and their target genes stored in the Ingenuity® 
Knowledge Base upstream regulator analysis was 
performed. Two statistical measures, standard in IPA, 
were used to detect potential TRs: an overlap P-value 
and an activation Z-score. First, the analysis examined 
how many known targets of each TR were present in 
our data set, resulting in an estimation of an overlap P- 
value. We set a threshold of an overlap P < 0.05 to identify 
significant upstream regulators. Second, the known effect 
(activation or suppression) of a TR on each target gene 
was compared with observed changes in gene expression. 
Based on concordance between them, an activation Z- 
score was assigned, showing whether a potential TR was 
in 'activated' (Z-score > 2), 'inhibited' (Z-score < —2) or 
uncertain state. 

Circos plots. We built Circos plots (34) to simultaneously 
visualize activated/inhibited upstream TRs (mainly con- 
sisting of TFs, Histone deacetylases (HDACs) and 
nuclear receptors), miRNAs and selected biological func- 
tions for each time point. Briefly, we first included known 
activated/inhibited upstream TRs based on SDE genes for 
each time point ('inferred TFs'). Next, we selected 
miRNAs targeting these inferred TFs from TarBase v. 6 
(http://www.microrna.gr/tarbase) (47). Then, we identified 
other non-TR target genes of these miRNAs and created a 
network representing inferred TRs, miRNAs targeting 
these TRs and mRNAs, targeted either by TRs and/or 
by selected miRNAs. Importantly, only mRNAs and 
miRNAs, which showed significant differential expression 
in at least one condition (compared with untreated 
controls), were included in the network. Next, the 
identified SDE mRNAs (that are targeted by TR and/or 
miRNA) were subjected to IPA for functional annotation. 
To increase readabihty of the final Circos plots, we 
combined biologically related Ingenuity® functions into 
12 functional categories. Finally, we visualized identified 
networks and connections between molecules and func- 
tional categories using Circos. 

Identification of regulatory loops. To identify regulatory 
loops involving TRs, miRNAs and their common 
targets, we first targeted our approach to build networks 
of previously inferred TRs and their corresponding SDE 
mRNA targets for each time point. Secondly, SDE 
miRNAs were overlaid onto the networks. Next, connec- 
tions were associated with miRNAs that share the same 
targets as TRs and which were at the same time targeted 
by the TRs (a specific type of FFL in the form of 'TF 
miRNA-targets', edges from both regulators to target). 
Finally, a subset of FFLs from the networks was extracted 
based on biological significance. All the connections of 
these FFLs were based on experimentally vaHdated inter- 
actions referenced in the Ingenuity® Knowledge Base, and 



the mRNAs and miRNAs were selected based on their 
connectivity and their significance of differential expres- 
sion at each time point. 

RESULTS 

In the current study, we analysed and integrated matched 
time-series miRNA and mRNA microarray data, which 
were generated from melanoma cells after IFN-y stimula- 
tion as outHned in Supplementary Figure SI. The previ- 
ously described dynamic miRNA data sets (40) were 
combined with newly generated and matched mRNA ex- 
pression levels from the same samples for selected time 
points. Optimal sampling time points had been determined 
before by monitoring expression changes in response to 
IFN-y stimulation of a small number of miRNAs and 
genes [data not shown and (40)]. Given a fixed budget, 
we mostly opted for microarray duplicates rather than trip- 
licates and could therefore include more time points. 
However, to increase sensitivity and to control for batch 
effects, we performed an additional third replicate for 24, 
48 and 72 h. Two replicates of paired miRNA-mRNA 
microarray expression data were obtained for the control; 
for 3, 12, 24, 48 and 72 h of IFN-y-treated samples; and for 
the Jll-pre- treated control (JII Ctrl). We have shown before 
that JII prevented the phosphorylation and thus activation 
of STATl TF and its downstream actions, and therefore 
served as an additional control (40). 

Analysis pipeline for identification of differentially 
expressed mRNAs and miRNAs 

The outHne of computational analysis steps as applied to 
all data sets is depicted in Figure 1 . A multi-step approach 
was performed to identify SDE miRNAs and mRNAs 
over time. First, we checked the quality of the time-series 
mRNA data sets using Partek® Genomics Suite. No 
outlier microarrays or batch effects were detected (data 
not shown). If not stated otherwise, we included the 
third mRNA replicates in the analysis. As illustrated by 
the correlation study (Supplementary Figure S2), reprodu- 
cibiHty for mRNA expression data was evaluated by 
calculating coefficients of determination (r^) as previously 
done for miRNA expression data (40). Average coefficient 
of determination with 95% confidence was equal to 
0.979 ± 0.005 for both replicates and 0.903 ± 0.004 for 
non-replicates, showing the robustness of expression 
data sets. After pre-processing and quahty control, a fil- 
tering step was included (as outlined in 'Materials and 
Methods' section), which reduced the number of con- 
sidered features for both data sets, and improved the 
estimated FDR. In total, 16193 mRNAs (out of 33 279) 
and 160 miRNAs (out of 1100) passed the filtering and 
were considered for further analysis. To identify the best 
possible method for analysing these data sets, we first 
compared three methods available as R/Bioconductor 
packages: 'limma', 'betr' and 'timecourse'. Empirical 
Bayes methods and linear modelHng ('limma') performed 
better than other methods in terms of flexibility (as there is 
no need to have equal numbers of experimental replicates) 
and in terms of FDR on permutated and simulated data 
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(Supplementary Data). Therefore, differential expression 
analysis was performed using the limma package. For 
further analyses, a threshold FDR < 0.001 was applied, 
resulting in 65 miRNAs and 6056 mRNAs SDE over all 
time points. Using the same model and contrasts between 
expression levels in IFN-y-treated samples versus un- 
treated controls, we identified SDE mRNAs for each 
time point separately. 

Dynamic expression changes of miRNAs are delayed 
relative to mRNA expression changes 

To obtain a global view of the behaviour of filtered and 
paired miRNA-mRNA expression data, principal compo- 
nent analysis (PCA) was performed on each of the two 
data sets. The first principal component of two independ- 
ent PCAs was plotted (Figure 2A), showing transcriptome 
evolution over time, with the horizontal axis representing 
variability in miRNome and the vertical axis representing 
variabiHty in mRNAs. Remarkably, the principal compo- 
nents of both mRNA and miRNA data showed strong 
and reproducible time effects and accounted for 39% of 
data variabiHty for mRNA and 58% for miRNA. Owing 
to the properties of the principal component, this repre- 
sentation shows the main variability of mRNA and 
miRNA expression levels of all samples. Early time 
points (3 and 12h) and both controls (untreated and JII 
pre-treated) cluster together implying that the overall 
changes of the transcriptome were comparable with the 
average variability between replicates. The behaviour 



suggests that miRNA expression changes after IFN-y 
stimulation were delayed with respect to mRNA expres- 
sion changes, which were observed already at earlier time 
points. Interestingly, miRNA levels continued to change 
until 72 h, whereas mRNA levels were not altered signifi- 
cantly after 48 h. This suggests that mRNA expression 
levels adapt faster to the cytokine stimulus, possibly to 
initiate a rapid inflammatory response, which is then 
followed by a second transcriptional wave where 
miRNAs are involved in the regulatory cascade to fine- 
tune and adjust the system responses in the form of 
feedback regulators. Such miRNA-mediated feedback 
and FFLs have been described as common network 
motifs (48), and this observation was further supported 
by our 'limma' differential expression analysis. Several 
contrasts (two-class comparisons) were generated 
through comparisons between IFN-y-treaded samples 
and untreated controls. The resulting number of signifi- 
cant mRNAs and miRNAs (FDR < 0.001) is shown in 
Figure 2B, where biggest effects of cytokine stimulation 
on the transcriptome were scored between 12 and 24 h for 
mRNAs and between 48 and 72 h for miRNAs. At 24 h, 
~45% of all significantly expressed mRNAs were 
up-regulated, whereas expression changes of miRNAs 
(>50%) occurred later (at 48 h). 

Hierarchical clustering based on the expression profiles 
of significantly regulated mRNAs and miRNAs is 
shown in Figure 3. For better readabihty, we only show 
the top 100 SDE mRNAs and ah 65 SDE miRNAs. 
Further analyses were performed using all significantly 




Figure 2. Projection of dynamic changes of the transcriptome and the miRNome of melanoma cells after IFN-y stimulation. (A) PCA visualizes the 
evolution of both data sets over time, with the vertical axis corresponding to the first principal component of mRNA data and the horizontal axis 
showing the principal component of miRNA data. The percentage of variabiHty in the data sets represented by each axis is shown. The dots represent 
sample duphcates (complete arrays) for the indicated stimulation times. (B) The number of SDE mRNAs and miRNAs with FDR < 0.001 for each 
condition compared with the untreated control is shown. The maximum gradient was observed between 12h/ctrl and 24h/ctrl for mRNAs and 
between 24h/ctrl and 48h/ctrl for miRNAs. 
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Figure 3. Heat maps representing the time evolution of the top 100 mRNAs and all 65 differentially regulated miRNAs detected through multi-class 
'limma' analysis (FDR < 0.001). Standardized expression values for each feature were reordered by hierarchical clustering, resulting in three 
pronounced clusters depicted on the right of each heat map. Each cluster contains member gene or miRNA profiles (grey fines) and mean expression 
values (dots). (A) Tfie top 100 significantly expressed and annotated mRNAs fall into three main groups (cluster A, B, C). Names of TFs involved in 
IFN-y signalfing are marked in red. Clustering was supported by bootstrapping using the 'clusterCons' package of R. Altered cluster annotations 
were only observed for two genes (A2M and APOBEC3G), which after bootstrapping were more fikely to belong to cluster C rather than B. (B) The 
majority of all SDE miRNAs belong to two clusters (a and c), with delayed up- or down-regulation after 24 h, respectively. Cluster b contains only 
two miRNAs (miR-125b* and miR-21*). 
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regulated mRNAs (FDR < 0.001) (Supplementary Figure 
S3). Based on this approach, we have identified three main 
clusters for mRNAs (Figure 3A), supported by a boot- 
strapping analysis using 'clusterCons' package of R (1000 
iterations of hierarchical clustering with Euclidian 
distance). Cluster A showed delayed mRNA up-regulated 
expression reaching a plateau at 48 h, whereas cluster B 
contained genes responding rapidly to IFN-y stimulation. 
The average profile showed up-regulation after IFN-y 
stimulation, followed by a decrease in expression down to 
almost baseline levels after 72 h. This latter cluster included 
well-known TFs involved in IFN-y-stimulated signalling: 
STATl and IRFl and also chemokines (CXCLIO and 
CXCLl 1), which are typically up-regulated by IFN-y treat- 
ment (49). Finally, cluster C included down-regulated 
mRNAs (after 24 h) that remained at low expression 
levels until 72 h. Figure 3B depicts miRNA expression 
changes for the selected time points. We have previously 
observed that miRNA expression patterns significantly 
change only after 48 h (40), and were able to further 
confirm these alterations using the herein described differ- 
ent analysis pipeHne. The JII control samples (72 h) 
demonstrate that mRNA and miRNA regulations were 
diminished by blocking the Jak/STAT signalHng cascade, 
indicating that the observed expression changes were 
caused by IFN-y-activated downstream regulators. It is 
known that IFN-y treatment activates a variety of target 
genes and among those several TFs (49). Our data again 
suggest that these activated TFs then start a second wave of 
transcriptional activations, which seems to include many 
miRNAs the expression levels of which subsequently 
increase with a certain time delay. 

Correlation and validation of dynamically regulated 
mRNAs and miRNAs 

Next, we identified positively and negatively correlated 
miRNA-mRNA pairs using CoExpress, an in-house bio- 
informatics tool (http://www.bioinformatics.lu/ 
CoExpress), which combines correlation analysis with 
miRNA target gene predictions (see Supplementary Data 
and Supplementary Table S3). A negative dynamic correl- 
ation between a given miRNA and a predicted target 
mRNA could be the first indicator of a potential direct 
interaction. Supplementary Figure S4A summarizes 
the number of co-expressed mRNAs for each of the 
co-expressed miRNAs (|r| > 0.95). The first nine miRNAs 
with most negatively correlated events (miR-31, -1308, - 
424*, -29b-l*, -23a*, -92a-l*, -29a, -22, -27a*) belonging 
to cluster a (Figure 3B) show a clear tendency to have more 
negatively correlated expressions with mRNAs, whereas 
miR-27b (cluster c) had the highest proportion of positively 
correlated events. Interestingly, increasing the stringency of 
the correlation threshold resulted in a higher proportion of 
negatively correlated events (Supplementary Figure S4B). 
Microarray-measured expression levels were vaHdated by 
qPCR for nine miRNA-mRNA pairs that were selected 
based on correlation, confirming expression patterns 
for all tested pairs (Supplementary Figure S5). We mainly 
focussed on two miRNAs, miR-29b-l* and miR-424*, 
that we have previously shown to be among the top 



10 dynamically up-regulated miRNAs after IFN-y 
treatment (40). 

CoExpress was used to identify negatively correlated 
and biologically interesting mRNAs that are known to 
be involved in tumourigenesis of melanoma or other ma- 
Hgnant tumours (50-53). TIAMl (T-cell lymphoma 
invasion and metastasis 1) and IGFBP5 (insulin-Hke 
growth factor-binding protein 5) were dynamically anti- 
correlated with miR-29b-l*. MMPl (matrix metallo- 
peptidase 1) and SOX5 [SRY (sex-determining region 
Y)-box 5] showed inverse correlations over time with 
miR-424*, while VCAN (versican) was negatively 
correlated with both miRNAs. Expression levels of all 
Jll-pre-treated samples remained largely unchanged over 
time for SOX5 and HDAC9, indicating that the observed 
changes in the miRNA levels were likely caused by the 
IFN-y-induced Jak/STAT signalling. In contrast, for 
MMPl and IGFBP5, mRNA levels in the Jll-treated 
control samples were also down-regulated, suggesting 
additional IFN-y-independent effects. Taken together, 
quahty measures as well as vaHdation of microarray 
results confirmed high quahty and reproducibihty of 
both data sets. In the following steps, we combined and 
analysed these time-series-derived data sets representing 
dynamic expression changes of the melanoma miRNome 
and transcriptome to identify interactions that only 
become visible over time. 

Functional annotation of dynamically regulated 
mRNAs and miRNAs 

To better understand which and how biological functions 
are affected by differentially regulated mRNAs and 
miRNAs over time, we performed functional annotations 
in Ingenuity® Systems (IPA) based on SDE mRNAs and 
miRNAs (each condition was compared with untreated 
controls). All 62 significant functional categories 
reported by IPA (adjusted P<0.01) are presented in 
Figure 4. To show the time evolution among the 
categories, the smallest adjusted P-values of member func- 
tions were taken for each category and each time point. 
Colouring was performed based on scaled logio-trans- 
formed adjusted P-values, with white colour correspond- 
ing to adjusted P>0.01 and dark grey corresponding to 
the smallest adjusted P-value for a given category among 
all time points. Visualization of the resulting dynamics 
allowed for comparison of functional enrichment 
between the five time points. We found five groups with 
temporal differences in enriched functional responses to 
IFN-y treatment. The first group included very early 
cellular reactions to IFN-y, with enrichment of functions 
mainly observed between 3 and 12h. In the second and 
largest group, functional enrichments peaked around 12 h 
after treatment and included several categories of immune 
and inflammatory responses. A third group showed 
smallest P-values for functions at intermediate time 
points 24 h after IFN-y treatment, while only four 
categories, including 'cancer' and 'post-translational 
modifications', had enrichment peaks at 48 h. We have 
previously observed that most miRNAs are being differ- 
entially expressed at this time point (40), suggesting a 
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Figure 4. Dynamic changes in inferred functional categories based on SDE mRNAs and miRNAs. The minimum adjusted P-values for member 
functions were combined to iUustrate dynamic changes in ah enriched functional categories obtained from IPA analysis ('n/s;: P> 0.05). The intensity 
of grey boxes represents scaled adjusted P-values (log-transformed) for each category: white — non-significant (>0.01), dark grey — smallest adjusted 
P-value for each category among the time points. 
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second wave of transcriptional changes, including the 
up-regulation of miRNAs after IFN-y treatment. 
Finally, several categories such as 'cellular movement' 
and 'cell cycle' were only enriched very late after 
cytokine stimulation. Taken together, functional annota- 
tions of dynamic miRNA and mRNA expression changes 
after IFN-y treatment showed that some categories were 
predominantly enriched at only one specific time point, 
emphasizing once more the importance of analysing 
time-series-derived data. 

Dynamic features of the interferon pathway 

To dissect the temporal behaviour of key players involved 
in interferon signalling in more detail, we analysed our 
data using the Ingenuity® Systems (IPA) program. As 
expected after stimulation with IFN-y, the interferon 
signalHng pathway was found to be highly significant. 



especially at early time points (P- values for 3h/ctrl and 
12h/ctrl were 4.8 x 10~^^ and 2.0 x 10~^^ respectively). 
Figure 5 illustrates the IFN-y sub-network as provided 
by the IPA Knowledge Base to which we connected ex- 
pressed miRNAs that have been described to target the 
pathway-associated mRNAs (TarBase v. 5 and TargetScan 
with total context score < —0.4) and were SDE in at least 
one time point. Most of the genes (65%) involved in the 
Jak/STAT signalling pathway were differentially ex- 
pressed in at least one time point, which is depicted by 
small adjacent bar charts showing the experimentally 
measured expression levels over five time points (expres- 
sion levels of the Jll-treated control samples are shown on 
the far right side). 

Although STATl needs to be phosphorylated to become 
activated, the gene itself is known to be up-regulated by 
IFN-y stimulation (54), which was also confirmed here. 
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Figure 5. Representation of the top canonical pathway 'interferon signalHng' detected when simultaneously analysing the mRNA and miRNA data 
sets with IPA. Parts related to IFN-a/P were removed, and SDE miRNAs targeting the detected genes of the pathway were added. Connections 
between miRNAs and their targets were estabhshed by IPA. Genes that were differentially expressed in at least one condition are marked with filled 
grey symbols. Expression changes at 3, 12, 24, 48 and 72 h with respect to untreated controls are shown as bar charts close to each molecule. For 
non-significant conditions, a fine is shown instead of a bar. The last bar on the far right always corresponds to Jll-treated control (72 h). Connections 
between main players of the signalHng pathway are depicted as Hues: relationships between miRNAs and target mRNAs are shown as thin Hues, 
whereas relationships between TFs and target mRNAs are presented by thicker lines either indicating activation or repression. 



2826 Nucleic Acids Research, 2013, Vol. 41, No. 5 



Three downstream targets of STATl (ARFl: ADP 
ribosylation factor 1, IRF9: interferon regulatory factor 9 
and TAPl: antigen peptide transporter 1) show an imme- 
diate up-regulation after 3h of interferon treatment, 
whereas most targets reached peak expression levels at 
24 h (IRFl: interferon regulatory factor 1, IRF9, TAPl, 
IFITMl: interferon-inducible transmembrane protein 1, 
IFI35: interferon-induced 35-kD protein and PSMB8: pro- 
teasome (3 subtype 8). Consequently, the transcriptional 
targets of IRFl (BCL-2 and BAX) became up-regulated 
in a second cascade, showing increased expression 
levels only at later time points. In Hne with previous obser- 
vations [reviewed in (49)], at 12 h, we observed a moderate 
up-regulation of the suppressor of cytokine signalling 
(SOCSl), a negative feedback regulator, which reduces 
the activation of the Jak/STAT pathway. The temporal 
expression profiles of miRNAs connected to this pathway 
reveal inverse correlations with expression levels of their 
possible target genes. For example, miR-301a-3p was 
down-regulated at late time points, with its target gene 
IRFl showing increased levels; miRNAs known to target 
BCL-2 were all down-regulated over time, whereas BCL-2 
was up-regulated at 72 h. BAK, a target gene of IRFl 
TF was not differentially regulated over time (white 
symbol) and miRNAs targeting this mRNA were either 
up-regulated (miR-29b-3p) or down-regulated (miR-27b- 
3p and miR-125b-5p). Interestingly, several mRNAs in 
the network showed similar expression profiles at 72 h in 
the IFN-y-treated cells and in the Jll-pre-treated cells 
(STATl, IRF9, IFI35 and IFITMl), suggesting that 
feedback inhibitory processes were progressively developed 
in IFN-y-treated cells recapitulating the phenotype 
observed in cells inhibited with JII. To summarize, 
visualizing dynamic expression data provides additional 
information about the interplay between miRNAs and 
mRNAs within the well-known IFN-y signalHng cascade. 

Integration of mRNA, upstream regulators and 
miRNA data 

Although initially designed and used for displaying com- 
parative genomic data (34), Circos plots have also been 
adapted to analyse mutations in cancer (55,56), meta- 
genomic data (57) and dynamics of TF regulatory 
networks (58). Here, we have applied Circos plots to inte- 
grate data sets from three different sources, and to 
our knowledge, this is the first time that data from 
miRNome and transcriptome were simultaneously 
combined with annotated functions (Figure 6). We chose 
to work with biological functions because it provides more 
robust results than working with individual genes. To 
decrease complexity of the graphs and to allow for 
better readability, we only show those TRs, miRNAs 
and inferred functions that were selected through the 
pipeHne described in 'Materials and Methods' section. 
Based on this integrative approach, we observed that 
miRNAs only become connected to TRs and biological 
functions at very late time points (48 and 72 h), indicating 
that they are activated by TRs downstream of the 
activated Jak/STAT signalHng cascade. Similar to what 
was seen in Figure 5, Circos plots representing 3h of 



IFN-y treatment and 72 h of Jll-control treatments 
could almost be superimposed, suggesting that new con- 
nections scored at later time points were indeed brought 
about by Jak/STAT signalHng. Although STATl was 
activated (by phosphorylation) after 15min of IFN-y 
treatment (40), it only became connected to annotated 
functions at the 12-h time point, while no direct inter- 
actions between STATl and selected miRNAs were 
detected at any time point. Furthermore, after 12h of 
IFN-y signalling, other TRs such as IRFl, nuclear 
factor of kappa Hght polypeptide gene enhancer in B- 
ceHs (NFkB) and Enhancer of zeste homolog 2 
(Drosophila) (EZH2) were predicted to be active based 
on expression changes of their regulated genes. IRFl, a 
TF activated by STATl, was actively 'communicating' 
throughout the experiment and could be involved in the 
downstream and delayed activation of several depicted 
miRNAs. At 48 and 72 h, we detected six connections 
between miRNAs and TRs. Additionally, the overall 
number of interactions between all players and inferred 
functions dramatically increased at these late time 
points. Again, this was not seen in the JII- treated 
negative control samples, implying that the estabHshed 
connections cannot be attributed to noise or mere 
'ageing' of cells in culture. Of note, miR-93-5p, one of 
three miRNAs depicted in the Circos plots, was not 
present in the heat map (Figure 3B), as the hierarchical 
clustering was generated based on multi-class 'limma' 
analysis and the adjusted P-value for miR-93-5p was 
just above the selected threshold (with FDR = 0.0013). 
To generate Circos plots, we used contrasts between the 
72-h time point and untreated samples in which 
miR-93-5p was detected as SDE, with an adjusted P- 
value of 0.0004 by two-class 'limma' analysis. 

Common transcriptional regulatory network motifs 
consisting of FFLs typically involve TFs, miRNAs and 
joint targets. Surveying possible instances of miRNA- 
mRNA-TF loops in our integrated time-lagged data sets 
aHowed for identification of several FFLs that only 
became fully active at specific time points (Figure 7). 
Five temporally interconnected FFLs appeared to be of 
special interest because they could control the expression 
of genes involved in biological functions that are relevant 
to our model (ceH adhesion, apoptosis and/or immune 
response or cell cycle). At 24 h, expression of ICAM-1 
(intercellular adhesion molecule- 1) gene appeared to 
be fine-tuned by a complex interplay between miR-21 
and RelA itself under the control of miR-193b. 
Simultaneously, the expression of the BCL-2 protein 
family member MCLl (myeloid ceH leukemia sequence 
1) gene was controlled by miR-29 and two TFs, the 
RelA/NF-KB complex and TP53. At 48 h, although the 
RcIA/NF-kB TF complex was no longer active, TP53 
still seemed to be involved in the miR-29-MCLl loop, 
but also controHed the expression of let-7a-5p miRNA 
and BCL2L1 (BCL-2 like-1) gene, another member of 
the BCL-2 protein family. Finally, at 72 h, an FFL 
involving let-7a-5p miRNA and KDM5B [lysine 
(K)-specific demethylase 5B] TR appeared to modulate 
the expression of CCNDl (cyclin Dl) gene. No active 
FFL were detected at 3, 12 and 72 h or in control 
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Figure 6. Circos diagrams showing dynamical dependency of the transcriptome, represented in the form of top biological categories (purple) and of 
inferred upstream regulators: transcription regulators (TRs: brown) and miRNAs (blue). Time points for 3, 12, 24, 48 and 72 h and for the 72-h 
Jll-treated control were compared with the untreated control, and the SDE molecules were analysed by 'limma', with contrasts (FDR < 0.001) as 
described in 'Materials and Methods' section. The width of the categories is related to the number of member mRNAs. A connection between a TR 
and a functional category means that this TF was detected as an activated or inhibited TR by upstream regulator analysis, and that its target genes 
grouping in the respective functional category were differentially expressed in at least one time point. A connection between a miRNA and a TR 
implies that this miRNA was SDE and was predicted to target the TR genes. Finally, a connection between a miRNA and a functional category 
indicates that one or several of its target genes of a differentially expressed miRNA belonged to the assigned category. The thickness of connecting 
lines illustrates higher number of targets of TRs or miRNAs within this functional category. 
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Figure 7. Graphical representation of regulatory sub-networks (extracted from IP A) includes three activated TRs (dark grey boxes), four genes 
(rounded boxes) and four miRNAs (ellipses). Activation time for each part of the sub-network is shown by arrows on top. 'mir' represents the 
immature form of miRNAs, whereas 'miR' denotes the mature form. Connections between molecules are presented with respect to experimental 
observations and IPA predictions. Black: molecules have correlated (anti-correlated in case of inhibition) expression profiles. Dotted arrows: mRNA 
profile of a target molecule is in concordance with the predicted activation state of the TF. Grey arrows: direct interaction was not observed, 
suggesting presence of cumulative effect of other regulators. 



samples treated with JII. Our data suggest that response to 
IFN-y in our cell model involves a discrete set of TRs that 
act sequentially and/or synergistically on cell stimulation. 

Taken together, the integrative analysis of data sets rep- 
resenting high-resolution time-series-derived mRNA and 
miRNA microarray data allowed us to detect dynamic 
interactions that would have been missed if only single 
time points were considered. Our biostatistical and bio- 
informatics analysis pipeline identified a very interesting 
time delay of the miRNome with respect to the transcrip- 
tome, and we speculate that most miRNAs become 
up-regulated only after the first round of transcriptional 
activation is completed. Furthermore, the incorporation 
of inferred biological functions added another level of 
complexity to this study, allowing for visualization of 
dynamic changes in functional systems. 



DISCUSSION 

Owing to their abihty to post-transcriptionally regulate 
gene expression of almost all genes, miRNAs are known 
to influence many cellular activities in healthy and 
diseased states. Because they are involved in key cellular 
processes, it remains essential to decipher more miRNA 
target genes and to examine how miRNAs are regulated, 
their temporal dynamic behaviour and their involvement 
in defined cellular functions. 

The current study was motivated by the question of how 
the functional interplay between mRNAs and miRNAs is 
regulated and changing dynamically. To explore the 
global temporal response to IFN-y treatment, we 
examined the expression levels of miRNAs and mRNAs 
in a time-series experiment using A375 melanoma cells. 
Time-series analysis, as opposed to comparison of 
multiple steady states, gives important insights into the 
causality of the observed interactions: while normally 
only connections between molecules are described, 
time-series data allow for addressing the direction of the 
interaction and its rate, and thus provide a better under- 
standing of cell dynamics (30,59). In a previous study (40), 



we performed a detailed investigation of the dynamic be- 
haviour of miRNAs over a wide time range after cytokine 
stimulation with IFN-y, which activates the TF STATl. A 
surprising finding of this former study was that all 
miRNA expression changes occurred with a delay only 
after 24 h. To find an explanation for this result and to 
identify dynamic regulatory networks, we performed a 
series of mRNA microarrays using the same RNA 
extracts. In addition to the identification of significantly 
regulated miRNAs and mRNAs over time, this approach 
allows for building and testing contrasts using the same 
Hnear model (45). Based on our benchmarking results of 
three methods ('betr', 'timecourse' and 'limma'), 'limma' 
proved to be superior in terms of FDR for permutated 
data sets and synthetic data. Using the 'limma' tool, we 
confirmed the previously reported 23 differentially 
regulated miRNAs (40), and further refined this data by 
detecting an additional 42 miRNAs with an FDR < 0.001 
(Figure 3B). 

Numbers of SDE genes and miRNAs (Figure 2B) and 
the heat maps (Figure 3) clearly revealed a delayed 
response of the miRNome to IFN-y stimulation with 
respect to the transcriptome. Interestingly, Pedersen 
et al. (60) have described two miRNAs that were 
modulated already after 30min after IFN-(3 stimulation 
of hepatocyte cells. Also, Kutty et al. (61) observed early 
regulation of miR- 155 after exposure to a tumour necrosis 
factor-a, interleukin- 1 (3 and IFN-y cytokine mixture. 
Here, we did not identify statistically significantly 
regulated mature miRNAs reacting that quickly to 
IFN-y in melanoma cells. In this context, we have 
recently confirmed that primary miRNA transcripts 
(pri-miRNAs) and precursor form (pre-miRNAs) (of the 
miR-29 family) are up-regulated well before the corres- 
ponding mature miRNAs after IFN-y stimulation of 
melanoma cells (62). The combined analysis of miRNA 
and mRNA data sets suggests that the IFN-y-initiated 
Jak/STAT signalling cascade transcriptionally activates 
other downstream TFs as well as other effector genes, 
which may in turn participate in up-regulation of 
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expression levels of a number of responding miRNAs in 
melanoma cells. 

On comparison of both data sets, we observed only 
three negatively regulated mRNAs at 3h, while 117, 
including STAT 1-3, IRFl-6, IFI16 and other IFN-A- 
related genes, were up-regulated. From 24 h onwards, 
more down-regulated mRNAs were scored, which chrono- 
logically was in good concordance with the elevated 
miRNA levels, indicating that there might be an inverse 
correlation for some of these miRNA-mRNA pairs. Three 
main dynamic profiles were detected among the top 100 
mRNAs (Figure 3A): late expression, early expression 
followed by repression and late repression. These profiles 
were also seen when analysing the expression profiles of 
all significant mRNA (Supplementary Figure S3). The 
65 SDE miRNAs, in contrast, presented only two major 
dynamic profiles (Figure 3B): delayed up-regulation and 
delayed down-regulation similar to our previous results, 
which included a more detailed analysis of temporal 
miRNA profiles (40). 

It is well accepted that identification of target mRNAs 
regulated by miRNAs is required to elucidate the exact 
role of individual miRNAs or groups of related 
miRNAs in a given cell. Several algorithms have been 
estabhshed for in silico predictions of target mRNAs [see 
review (25)]. However and as mentioned before, there is no 
efficient algorithm that reliably predicts all targets with a 
minimal number of false positives. A straight-forward 
approach to improve target gene predictions could be 
global correlation analyses of miRNAs and experimen- 
tally matched mRNA expression patterns in combination 
with standard target gene prediction algorithms at least 
for those interacting pairs where miRNAs cause a 
measureable decrease in mRNA levels. For this, we used 
the in-house-developed tool CoExpress to build a 
miRNA-mRNA correlation map, to compare negatively 
correlated miRNA-mRNA pairs with TargetScan predic- 
tions and to experimentally confirm interactions extracted 
from TarBase. Using this approach, we detected 398 
negatively correlated predicted targets for 21 miRNAs 
(Supplementary Table S3 based on TargetScan), and we 
were able to validate 14 selected miRNA-mRNA pairs by 
qPCR (Supplementary Figure S5). 

Overall, integrating time-series miRNA and mRNA 
data sets provides valuable information not only for iden- 
tification of possible miRNA target genes but also for the 
elucidation of dynamic changes of the underlying biolo- 
gical processes. Using Circos plots, we visualized specific 
interactions between upstream TRs, miRNAs and also 
mRNAs that were categorized in a set of biological func- 
tions (Figure 6). Circos facihtates the integration of 
different data sets, thus providing a more holistic view 
of the evolving processes. Ebert and Sharp have recently 
summarized evidence that suggests miRNAs to actively 
confer robustness to biological processes by dampening 
and/or increasing cellular responses to internal or 
external stimuli (24). This random noise in transcription 
rates, and as such in expression levels of genes and 
proteins, may in part be kept within certain boundaries 
by the action of miRNAs. Our data on temporal expres- 
sion changes within the miRNome and transcriptome as 



well as on TRs support this notion: in response to an 
external stimulus (IFN-y), many TRs and other mRNAs 
react rapidly with measurable expression level changes, 
whereas miRNAs react with a considerable time delay. 
This tentative 'second wave' of responses then brings 
up-regulated miRNAs into play, which mostly down- 
regulate their respective target genes, and thus control 
and reduce inordinate cellular reactions. 

Several interactions in the form of negative or positive 
FFLs have already been established between TFs and 
miRNAs (7,19,63). Very recently, Gerstein and collabor- 
ators described a comprehensive architecture of the 
human transcriptional regulatory meta-network derived 
from integrating data from the Encyclopedia of DNA 
Elements project with genomic and protein information 
(64). The authors found that this meta-network exhibited 
a high enrichment in FFLs, showing the importance of 
this motif in transcriptional regulatory process. Looking 
for specific FFLs, which include TFs regulating the ex- 
pression of both a miRNA and its target mRNA at the 
same time (FFL type: 'miR TF-targets', with edge 
from both regulators to target), we found several regula- 
tory relationships that could be of interest in our biolo- 
gical system. By combining time-series stimulation of gene 
expression by TFs with delayed transcriptional repression 
by miRNAs, FFL motifs may create a toggle-switch mech- 
anism by which the cell could timely coordinate responses 
to IFN-y stimulation. Chalancon et al. (20) recently stated 
that to gain a better understanding of gene regulatory 
events in response to environmental changes, it is neces- 
sary to measure expression changes at highest possible 
resolution and under many different conditions over 
wider time ranges. Although our study describes a rela- 
tively small number of matched data sets, it represents the 
first attempt to move to more quantitative time-resolved 
data and could be considered as a proof of principle for 
more extensive studies in future incorporating more con- 
ditions, time points and cell systems; the herein established 
computational analysis pipeline can easily be adapted 
to additional data sources or requirements. Knowledge 
from integrated data on individual miRNAs, famihes of 
miRNAs and eventually on the entire miRNome together 
with dynamic and matched transcriptome data will con- 
tribute to a more comprehensive view of biological 
systems, their regulation and their behaviour over time. 
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